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SECTION  I 


INTRODUCTION 


In  the  last  two  decades  there  has  been  remarkable  technical  progress 
in  the  fields  of  electronics,  in  general,  and  data  processing,  in  particular. 
In  the  same  period,  the  new  area  of  computational  fluid  dynamics,  a  branch 
of  numerical  analysis,  has  experienced  a  growth  comparable  to  that  of  com¬ 
puter  technology.  Among  the  many  applications  of  computational  fluid 
dynamics,  the  numerical  solution  of  the  Navier-Stokes  equations  has  chal¬ 
lenged  a  large  number  of  engineers  and  scientists,  due  to  the  capability  of 
these  partial  differential  equations  of  correctly  modeling  most  of  the  very 
interesting  phenomena  associated  with  fluid  viscosity  and  compressibility 
(e.g.  shocks,  shock-boundary-layer  interaction,  separation,  stall,  etc.). 

The  present  work  is  concerned  with  the  problem  of  low  speed  viscous 
flow  past  an  arbitrary  two-dimensional  body,  for  which  all  compressibility 
effects  are  negligible.  Even  for  the  case  of  the  incompressible  Navier- 
Stokes  equations,  the  number  of  numerical  techniques  and  their  applications, 
available  in  the  technical  literature,  is  so  high  that  no  systematic  survey 
will  be  provided  here.  For  the  present  purpose,  only  a  few  typical  examples 
will  be  mentioned:  the  pioneering  work  of  Burggraf^  concerning  the  now 

2_  A 

classical  driven  cavity  problem;  the  studies  of  Davis  and  his  co-workers 

for  laminar  flow  past  several  semi-infinite  bodies  at  moderate  to  high 

values  of  the  Reynolds  number;  the  analysis  of  Mehta  and  Lavan5  about  the 

starting  vortex,  separation  and  stall  of  a  lifting  airfoil.  However,  most 

2-5 

numerical  techniques  are  limited  to  particular  geometries,  for  which  a 
coordinate  transformation,  mapping  the  body  surface  into  a  coordinate 
line,  or  equivalently  the  potential  flow  solution  could  be  obtained 


analytically.  Two  approaches  appear  at  present  very  promising  for  removing 
such  a  difficulty  and  providing  viscous  flow  solutions  about  arbitrary  con¬ 
figurations:  The  Finite  Element  Method  and  the  numerical  generation  of 
body  oriented  coordinates.  In  particular,  Thompson  et  al^  '*'0  have  developed 
and  improved,  throughout  the  years,  a  numerical  technique  for  solving  the 
time  dependent  Navier-Stokes  equations  past  one  or  more  arbitrary  two 
dimensional  bodies:  First,  they  generate  an  appropriate  body-fitted  coor¬ 
dinate  system  which  maps  the  flow  domain  of  arbitrary  shape  in  the  physical 
plane  into  a  rectangle  in  the  transformed  plane.  Second,  they  solve  the 
unsteady  Navier-Stokes  equations  in  the  transformed  plane  (coordinates) 
by  means  of  an  implicit  time  marching  numerical  technique.  A  point  SOR 
(Successive  Over  Relaxation)  iterative  procedure  is  used  at  each  new  time 
level  in  order  to  solve  for  the  nonlinear  terms  and  the  elliptic  part  of 

the  equations,  explicitly.  This  approach  has  been  shown  to  be  applicable 

8  9  10 

to  both  the  vorticity-stream  function  and  pressure  velocity  ’  formulations 
of  the  time-dependent  Navier-Stokes  equations,  for  laminar  as  well  as  tur¬ 
bulent  flows.  Further,  it  has  been  proved  very  reliable  in  modeling  highly 
separated  flows  around  stalled  airfoils^.  However,  its  use  for  design 
purposes  is  severely  limited  by  its  computational  inefficiency.  In  part¬ 
icular,  when  it  is  used  to  provide  a  steady-state  flow  solution  by  following 
the  asymptotic  time  decay  of  an  unsteady  flow  phenomenon,  the  intrinsic 
inefficiency  of  point  iterative  methods  compounds  to  that  of  time  dependent 
approaches. 

The  aim  of  the  present  research  is  to  develop  an  efficient  numerical 
procedure  for  solving  the  steady-state  Navier-Stokes  equations  past  an 
arbitrary  two  dimensional  body,  by  combining  the  transformation  of  Thompson 
et  al6,7  with  a  numerical  technique  more  efficient  than  the  point  SOR  method. 


2 


Recently,  many  researchers  have  applied  a  number  of  ADI  (Alternating 

Direction  Implicit)  techniques  to  the  numerical  solution  of  the  Navier- 

2-4  15 

Stokes  Equations.  In  particular,  Davis  and  his  co-workers  *  ,  Briley 

and  McDonald^  ^  and  Beam  and  Warming^  have  obtained  considerable  success 

with  such  a  technique.  The  two  major  advantages  of  the  Linearized  Block 
13 

Implicit  methods  (like  the  ADI)  are  the  (quasi) linearization  of  the 

governing  equations,  which  eliminates  any  need  of  iterations  at  any  time 

level,  and  the  presence  of  only  block-tridiagonal  matrices,  whose  direct 

16 

inversion  is  performed  very  efficiently  by  block  Gaussian  elimination 

In  the  present  study  an  ADI  technique  will  be  used  to  solve  the  vorticity- 

stream  function  Navier-Stokes  equations  in  the  transformed  plane  after 

mapping  the  flow  field  around  an  arbitrary  airfoil  into  a  rectangle  by 

means  of  the  transformation  of  Thompson  et  al^’^.  Since  only  the  steady 

state  solution  is  of  present  concern,  the  stream  function  equation  is 

parabolized  by  adding  to  it  a  relaxation-like  time  derivative,  according  to 
2 

Davis  .  The  vorticity  equation  is  then  (quasi) linearized  and  the  two 

(equations)  are  solved  as  a  coupled  set  of  finite  difference  equations 

by  means  of  the  Douglas  and  Gunn‘S  ADI  technique.  The  incremental  approach 

12  14 

of  Briley  and  McDonald  also  used  by  Beam  and  Warming  and  by  Davis  and 
Hill'*'”’  is  used  at  the  second  sweep  of  the  ADI  procedure,  in  order  to  min¬ 
imize  computer  storage. 

The  coordinate  transformation^’ ^ ,  employed  here,  introduces  a  cut  in 
the  physical  plane,  which  is  mapped  into  the  two  vertical  sides  of  the 
integration  rectangle  in  the  transformed  plane.  Therefore,  the  additional 
difficulty  of  periodic  boundary  conditions  in  the  horizontal  direction  had 


1  O 

to  be  dealt  with.  To  this  end,  the  method  of  Ahlberg  et  al.  for  inverting 
a  tridiagonal  periodic  matrix  has  been  generalized  to  the  present  case  of 


3 


a  two-by-two  periodic  block  tridiagonal  matrix.  All  the  details  of  the 
algorithm  and  the  results  of  its  application  to  a  simple  model  problem 
are  given  in  the  Appendix. 

The  present  numerical  technique  has  been  applied  to  three  problems. 
First,  a  simple  Poiseuille  flow  has  been  computed  in  order  to  verify  the 
second  order  accuracy  of  the  method  versus  an  exact  analytical  solution. 
Second,  the  classical  driven  cavity  problem^  has  also  been  solved  to 
further  verify  the  proposed  algorithm  in  the  case  of  a  truly  two-dimensional 
flow  problem.  Finally,  the  flow  past  a  NACA  0012  airfoil  has  been  computed 
to  demonstrate  the  capability  of  simulating  the  viscous  steady  flow  past 
an  arbitrary  two-dimensional  body. 


SECTION  II 


] 

j 

GOVERNING  EQUATIONS  AND  COORDINATE  TRANSFORMATION 


The  governing  equations  are  the  nondimensional  vorticity  stream  function 

Navier-Stokes  Equations,  with  a  relaxation-like  time  derivative,  ,  added 

2 

to  the  stream  function  equation  in  order  to  parabolize  it: 
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Equations  (1)  and  (2)  constitute  a  set  of  parabolic  (in  time)  partial  dif¬ 
ferential  equations  which  can  be  solved  numerically  by  means  of  a  time  marching 
ADI  procedure.  However,  it  is  important  to  realize  that,  since  equations  (1) 
and  (2)  are  not  the  unsteady  Navier-Stokes  equations,  a  correct  description 
of  the  transient  is  not  provided  and  only  the  converged  solution  will  have 
physical  meaning. 

In  order  to  solve  equations  (1)  and  (2)  for  flow  past  an  arbitrary  two- 

dimensional  body  (e.g.  an  airfoil),  the  transformation  of  Thompson  et  al.^’^ 

is  used  to  generate  numerically  a  system  of  body  oriented  coordinates.  With 
6  7 

this  transformation  ’  the  flow  field  in  the  physical  plane,  comprised  between 
a  circle  of  radius  equal  to  ten  (chord  lengths)  and  the  airfoil,  is  mapped 
into  a  rectangle  in  the  transformed  (£,  n)  plane.  The  airfoil  and  the  circle 
are  mapped  into  the  lower  and  upper  sides  of  the  integration  rectangle 
respectively  and  the  two  sides  of  an  arbitrary  cut,  connecting  the  trailing 
edge  of  the  airfoil  to  the  outer  circle,  are  mapped  into  the  two  vertical 
sides  of  the  rectangle.  In  this  way,  the  two  vertical  boundaries  in  the 
transformed  plane  correspond  to  the  same  physical  line  and,  therefore,  periodic 


boundary  conditions  in  the  horizontal  (£)  direction  are  required.  The 
transformation  is  provided  by  a  set  of  two  elliptic  partial  differential 
equations  which  are  discretized  and  solved  numerically  by  means  of  a  point 
SOR  method^* The  step  sizes  in  the  transformed  plane  (£,  n)  are 
arbitrary,  since  they  cancel  out  in  the  coordinate  transformation  finite 
difference  equations,  and  are  both  taken  equal  to  one,  for  convenience**’2. 
Further  details  are  given  in  References  6  thru  10.  The  transformation  of 
Thompson  et  al.**’2  has  been  used  satisfactorily  in  several  numerical  solu¬ 
tions  of  viscous  and  potential  flows  in  regions  containing  any  number  of 

8—10 

arbitrary  two-dimensional  bodies  and  can  be  extended  to  three-dimensional 
configurations.  Its  two  major  limitations  are  due  to  the  approximation  intro 
duced  by  imposing  the  free-stream  boundary  conditions  at  a  finite  distance 
from  the  body  and  to  its  inability  of  removing  exactly  sharp  edge  singu¬ 
larities.  Whereas  the  boundary  condition  approximations  are  considered  suf¬ 
ficient  for  the  present  study,  the  second  limitation  has  been  removed  by 
considering  an  airfoil  with  a  rounded  trailing  edge. 

In  the  transformed  coordinates  the  governing  equations  (1)  and  (2) 
become : 
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where  J,  a,  6,  Y»  o.  T  are  the  Jacobian  and  the  scale  factors  of  the  coor¬ 
dinate  transformation,  see  Reference  10  for  their  analytical  expressions. 
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The  no-slip  and  zero  injection  boundary  conditions  at  the  surface  of 

Q 

the  airfoil  are  given  in  the  transformed  plane  as: 


>P  <C,  0)  =  0 

(5) 

*  (e,  o)  =  o 

(6) 

The  free  stream  conditions,  imposed  on  the  circle  enclosing  the  computational 
flow  field,  are: 


u)(C ,  nM)  =  0  (7) 

i(> C C »  r^)  =  yc  cosQ  -  xc  sin0  (8) 

where  0  is  the  angle  of  attack  of  the  free-stream  flow,  and  xc>  y are  the 
physical  coordinates  of  the  circle  corresponding  to  £  and  in  the  trans¬ 
formed  plane,  being  the  height  of  the  integration  rectangle,  equal  to  the 
number  of  gridpoints  (in  the  n  direction)  minus  one. 

Finally,  the  coordinate  transformation  introduces  the  following  addi¬ 
tional  (nonphysical)  periodic  boundary  conditions: 

(9) 

(10) 


iKSm,  n)  =  ^(0,  n) 


and 


u>UM,  n)  =  «(0,  n) 


where  is  the  width  of  the  integration  domain,  equal  to  the  number  of 
gridpoints  (in  the  C  direction)  minus  one.  As  previously  mentioned, 
boundary  conditions  (9)  and  (10)  produce  periodic  two-by-two  block  tri¬ 
diagonal  systems  in  the  second  sweep  of  the  ADI  solution  procedure.  Such 

a  difficulty  has  been  resolved  in  the  present  study,  where  the  Algorithm 
18 

of  Ahlberg  et  al.  for  solving  periodic  tridiagonal  system  has  been 
generalized  to  the  use  of  periodic  systems  of  two  coupled  tridiagonal  equa¬ 
tions  (see  the  Appendix) . 
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SECTION  III 


NUMERICAL  METHOD 


Equations  (3)  and  (4)  are  expressed  in  finite  difference  form  and  solved 

18 

numerically  by  means  of  the  Douglas  and  Gunn  ADI  procedure  as  follows: 

First,  the  nonlinear  convective  terms  in  the  vorticity  equation  are  (quasi)- 
linearized  and  the  time  derivative  are  replaced  by  finite  differences  to  give: 

,  n+1  n,  .  /,n+l  n  .  ,  n  n+1  n  nw  T 

(o)  -  to  )/At  +  oi_  +  ijt^  ui^  - 

.n+1  n  n  n+1  n  n,  , ,  .  n+1  Oo  n 

•(*5  un  “n  “  “n)/J  "  (“  “«  "  20  % 


n+1  n+1  n+1 ,  . 

+  >  ui  +  o  w  +  t cj _  )/Re  *  0 
nn  n  S 


n+1  ,n  ,  ,n+l  ,  ,n+l  ,  ,n+l  ,  n+1  /  n+1  ,n. 

a  \p  -  26  +  y  d'nn  +  o  >i'n  +  +  o>  “  ~  i>  ) / At  (12) 

2 

where  the  J  dividing  a,  8,  Y»  cr  and  i  has  been  omitted  for  convenience. 

Note  that  in  equations  (11)  and  (12)  all  the  linear  terms  are  expressed  at 
the  new  time  level  tn+1  *  tn  +  At  implicitly  except  for  the  mixed  derivatives 
which  are  expressed  (explicitly)  at  the  old  time  level  tn,  see  References 
14  and  15.  Alternatively,  equations  (3)  and  (4)  can  be  linearized  in  time 
according  to  a  Crank-Nicolson  averaging  to  give: 

,  n+1  nw..  ,  /.n+1  n  n  n+1.,,,,  ,,n+l  n  .  ,n  n+1. 

n  s,  n  £  t,  n  £  n 


ra  ,  n+1  n  ,  n  y  .  n+1  n  . 

'IT  Kr  +  “ft  *  26  id  +  T  (oi  +ai  ) 

2  ££  CC  Cn  2  nn  nn 


+  £  („.n+l  x  n.  t  /  n+1  n,,  n 

T  1  +  0)  )  +  7T  (o)  +  0).)]  “  0 


and 


a  /,n+1  .  ,n  .  ,n  .  y  ,  ,n+l  ,n  .  a  .,,n+l  ,n, 

7  (*rr  +  *rr)  ~  20  *r„  +  *  +  ♦„.>  +  7  O*  +  O 


n+1 


n+1 


K££ 


r££'  £n  2  'Tnn 


,  t  / . n+ \  . n  v  , 

+  + 


nn 

n+1 


a> 


+  0) 


nn 

i n+1  i 

t _ -...». 

At 


(14) 


which,  except  for  the  explicit  mixed  derivatives,  are  second  order  accurate 

in  time.  The  two  step  ADI  procedure  of  Douglas  and  Gunn^2  is  then  applied 

to  equations  (11)  and  (12)  (or  equivalently  to  equations  (13)  and  (14)). 

In  the  first  sweep,  the  solution  (indicated  by  a  *)  is  advanced  in  time 

by  evaluating  implicitly  only  the  n  derivatives  and  the  source  term  (wn+1) 

in  the  stream  function  equation,  and  evaluating  all  the  other  terms  at  the 

old  time  level  tn.  Equations  (11)  and  (12)  thus  become: 

* 


, n  n 
At  “  (~T  Re  Un  "  unn  +  “j 


<P 


n 


UJ  .  ,  n 

aT  +  (au)££ 


20io"  +  Tw”)/Re 

£n  £ 


and 


*  ip  *  * 

w  ~  77  +  +  y\p 

At  n  nn 


fe  •  T^£  “  “*££  +  28^n 


(15) 


(16) 


Equation  (15)  and  (16),  after  that  all  the  derivatives  are  replaced  with 
(second-order-accurate)  central  finite  differences  give  a  coupled  set  of 

linear  tridiagonal  equations  of  the  type: 

*  *  *  *  *  * 
al 

and 


’j-l 

+  blj 

“j 

+  ClJ 

VI 

+  dlJ 

Vi 

+  elj 

♦j 

+  flj 

Vi  '  hlj 

(17) 

* 

+  b2j 

* 

* 

* 

* 

* 

j-l 

"j 

+  c2j 

Vi 

+  d2J 

♦j-l 

+  e2j 

♦j 

+  f2j 

♦j+1  “  H2j 

(18) 

In  equations  (17)  and  (18)  the  subscript  i  (indicating  the  longitudinal  £ 
location)  has  been  dropped  for  convenience,  j  varies  from  2  to  J-l  and  all 
the  coefficients  are  known  and  can  be  obtained  straightforwardly  from 
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equations  (15)  and  (16).  At  all  locations  (i)  equations  (17)  and  (18)  are 
solved  very  efficiently  by  Gauss  block-tridiagonal  reduction^  thus  pro- 

it  it 

viding  all  the  id.  .  and  ip.  .  values  (i  ■  and  j  ■  1,...J).  Details 

1»J  J-jJ 

of  the  boundary  conditions  will  be  provided  later.  In  the  second  sweep 
of  the  ADI  procedure,  the  final  ('J»n+'*',  u)1^)  solution  is  obtained  by 
evaluating  the  £  derivative  and  the  source  term  wn+2  (in  equations  11  and 
12)  implicitly  and  the  n  derivatives  explicitly  from  the  first  sweep  * 
solution,  that  is: 

u>n  .  n+l  a  n+1  %  ,n+l 

—  +  (~J  “  R^)a)£  "  Ii“« 


-  -  — r  (i p  -  ip  )  +  — r  (a)  -  ui  )  +  ■—  +  (yw  +  Ou) 

J  n  n  J  n  n  At  1  nn  n 


-  2(3  gj“  )/Re 


u)n+1  +  t*"+1  +  a*"*1  -  *n+1/At  -  -  o»*  -  yp*  +  26*”  -  *n/At  (20) 

£  ££  n  nn  £n 
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According  to  Briley  and  McDonald  ,  equations  (19)  and  (20)  are  replaced 
by  the  following  ones  obtained  by  subtracting  equation  (15)  from  (19)  and 

(16)  from  (20),  that  is: 

in  n  . 

-ip  u  *  n 

w»/nTv“  a—  n  —  u)  —  gj  .  „ 

At  +  ^  J  "  Re^  U£  Re  U££  J  ^£  "  At 


w  +  t*  +  a*  - 

At 


'Af‘  V 


where  the  new  variables 
n+1  n 

w  =  (l)  -  u> 


7  . n+1  , n 

^  ^  —  ifi 


have  been  introduced  for  convenience. 


Equations  (21)  and  (22),  after  that  all  the  derivatives  are  replaced 
by  central  finite  differences,  become  a  set  of  coupled  linear  tridiagonal 
equations  of  the  type 

ali  w±_1  +  blt  wi  +  cl±  wi+1  +  dl±  i1_1  +  eli  +  fl±  wi+1  «  hli  (25) 


a21  toi_1  +  b2t  oi1  +  c21  to±+1  +  d2t  ji1_1  +  e2±  ^  +  f2t  i|»±+1  »  h2it  (26) 

where  the  subscript  j  is  now  dropped  for  convenience  and  all  the  al^  thru 

h2^  coefficients  are  known.  Equations  (25)  and  (26)  constitute,  for 

i  =  1,...I,  a  system  of  21  coupled  tridiagonal  equations  subject  to 

periodic  boundary  conditions,  so  that  for  i  *  1,  and  ijl  ^  are  replaced 

by  uij  and  ^  and  for  i  =  I,  and  are  replaced  by  and 

This  system  is  solved  very  efficiently  by  means  of  the  algorithm  presented 

in  the  Appendix  for  all  rows,  i.e.,  for  j  =  2,...J  -  1.  It  is  worth  noting 
*  * 

that  the  u).  and  ^  .  values  appearing  in  the  coefficients  hi  and  h2  are 

1 » J  1  *  J  J  j 

not  needed  for  the  evaluation  of  w,  and  ij),  ,  in  any  successive  row, 

*■*  J  >  J 

Am  A 

Therefore,  the  same  arrays  are  used  to  store  w.  and  u>,  .  and  ^  and 

1*3  *  >  J  J 

The  solution  at  the  new  t  time  can  now  be  evaluated  as: 


n+1  n  .  - 

“i,J  “  “i.j  Wi J 


(27a) 


,n+l  ,n  .  7 

*i.j  ‘  *i.J  +  *i,J 


(27b) 
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for  i  ■  !,...!  and  j  ■  2,...J  -  1  and  as 


n+1  * 

“i.j  “  “i.J 

(28a) 

and 

n+1  * 

♦i.j  “  *i.j 

(28b) 

for  i  ■  1, . 

. . I  and  j  -  1  or  j  -  J. 

The  whole  process  Is  then  repeated  until  convergence. 


A.  Boundary  Conditions  for  the  *  Solution 

In  the  first  sweep  of  the  numerical  procedure  equations  (17)  and  (18) 
have  to  be  solved  at  every  longitudinal  location,  subject  to  boundary  con¬ 
ditions  (5-8). 

Equations  (5),  (7)  and  (8)  are  immediately  imposed  as 

-  0  (29) 

**  -  0  (30) 

w*  -  0  (31) 


Equation  (6)  can  be  satisfied  in  several  different  ways.  The  following 

five  approaches  were  used  in  the  present  study.  Three,  four  and  five-point 

one-sided  finite  difference  representations  of  equation  (6)  give  respectively 

* 


(the  step  size  is  equal  to  one  and 


4 


(9 


2 

0* 


18^*)/2 


1^6  tjj* 

5*3  4 


12*;  +  ib*; 


is  eliminated  due  to  equation  (29)): 

(32) 

(33) 

(34) 
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■k. 


LI 


For  the  case  of  Cartesian  coordinates  (n  =  y,  y  =  1,  a  =  0)  a  linear  shear 
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flow  was  also  assumed  near  the  body  surface  which  gives  the  following 
equation: 

*  1  *  * 

a>l  +  j  a>2  +  3^2  **  0  (35) 

* 

Finally,  a  central  finite  difference  for  ^  was  used,  consistently  with 
the  overall  numerical  procedure,  i.e., 

^2  “  *0  “  0  (36) 

* 

However,  the  value  of  ^1^,  interior  to  the  body  was  not  available  and  had 
to  be  evaluated  somehow.  To  this  purpose  the  steady  state  stream  function 
equation  along  the  n  =  0  (j  =  1)  line,  where  all  the  £  derivatives  vanish 
identically,  is  easily  seen  to  provide 


*  * 

yip  +o)  =  0  , 

Tin 

which,  in  finite  difference  form,  becomes 

Y1(’J»2  "  2<|/*  +  i|Jq)  +  w*  -  0 

This,  combined  with  equations  (29)  and  (36)  finally  gives 
*  *  * 


(37) 


(38) 


2^2  +  u>1 


0 


(39) 


where  y^  is  evaluated  by  a  three-point  extrapolation  from  y and  y^ . 
Any  of  equations  (32)  thru  (35)  or  (39)  is  easily  satisfied  in  the  block 
tridiagonal  inversion  of  equations  (17)  and  (18).  Such  an  inversion  is 
similar  to,  and  simpler  than,  that  given  in  the  Appendix.  In  particular, 
it  provides  recursion  relations  of  the  type 
*  * 


Rl. 

J 

wj-l 

+  s1j 

♦j-1 

+  T1 

R2 . 

* 

U)  .  i 

+  S2 , 

* 

+  T2 

j 

j-1 

j 

i 

j-1 

(40) 

(41) 
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fv  - 

i  v  L-  eii , 


where  Rl^  thru  T2^  are  given  in  terms  of  Rl^+^  thru  T2j+1  and  can  be  easily 
determined  for  j  ■  J  -  1,  J  -  2, . . .  ,2  since,  from  boundary  conditions  (30) 
and  (31): 

R1J  =  *  slj  =  T1j  =  R2j  =  S2j  =  T2j  =  0  (42) 

Any  of  equations  (32)  thru  (35)  or  (39) ,  together  with  (40)  and  (41)  easily 

*  * 
provides  a  relation  for  in  terms  of  known  coefficients,  and  all  and 

* 

ip  j  can  finally  be  evaluated  by  means  of  the  recursion  formulas  (40)  and  (41) . 


SECTION  IV 


RESULTS 

The  present  numerical  technique  was  applied  to  three  different  pro¬ 
blems:  a  simple  Poiseuille  Flow,  for  which  the  second  otder  .accuracy  of 
the  method  could  be  verified  versus  the  exact  analytical  solution;  the 
driven  cavity  problem  for  which  a  fully  two-dimensional  solution  could  be 
compared  with  available  numerical  results  and  finally  flow  past  a  NACA  0012 
airfoil  in  order  to  assess  the  capability  of  the  present  method  to  achieve 
its  goal  of  computing  arbitrary  two-dimensional  flow  fields. 

A.  Poiseuille  Flow 

For  laminar  steady  flow  inside  a  two  dimensional  channel,  at  any  loca¬ 
tion  x,  the  longitudinal  velocity  profile  is  parabolic  and  the  normal  velo¬ 
city  is  zero  everywhere.  For  convenience,  the  maximum  velocity  at  the 
center  of  the  channel  (y  =  h)  was  taken  to  be  equal  to  h.  The  exact  analy¬ 
tical  solutions  for  the  vorticity  and  the  stream  function  in  the  lower  half 
of  the  channel  (0  <_  y  h)  are  therefore  given  as: 


uj  =  -2  +  2^ 
n 

(A3) 

and 

0  , 3 

* =  y  ■ 

(AA) 

Two  numerical  solutions  were  obtained  by  using  5  gridpoints  in  the  longitu¬ 
dinal  x  (x  =  £)  direction  and  11  and  21  points  in  the  normal  y  (y  =  n) 
direction,  so  that  h  was  equal  to  10  and  20  respectively  (Ay  H  1).  The 
one  dimensional  nature  of  the  solution  was  captured  perfectly  thanks  to  the 
periodic  boundary  conditions  in  the  £  direction.  The  vorticity  at  the  wall 
was  found  to  be  equal  to  -  1.9900  and  -1.9975  (for  11  and  21  gridpoints 
respectively),  when  using  the  point  image  boundary  condition  of  Burggraf1 


equation  (39),  (y^  =  1»  f°r  this  case)  thus  verifying  the  second  order 
accuracy  of  the  method.  The  other  four  boundary  conditions  for  <|»  =  ^  =  0 

at  the  wall  were  also  used  and  found  totally  satisfactory.  Actually, 
equations  (33),  (34)  and  (35)  all  replicated  the  exact  (analytical)  results. 
It  was  mentioned  that  an  "implicit"  and  a  Crank  Nicolson  time  splitting  of 
the  governing  equations  could  be  used,  see  equations  (11)  thru  (14);  both 
approaches  have  been  implemented,  found  to  be  unconditionally  stable  for 
this  model  problem  and  have  produced,  at  convergence,  identical  results. 

B.  Driven  Cavity  Problem 

The  present  algorithm  was  also  applied  to  the  classical  square  cavity 
problem'*'.  For  this  case  the  boundary  conditions  are  nonperiodic  in  both  x 
and  y  directions.  The  stream  function  tp  is  prescribed  to  be  zero  on  all 
walls  of  the  unit  square  flow  field  and  the  derivative  of  ip  in  the  direction 
normal  to  the  wall  is  equal  to  one  on  the  top  side  of  the  square  and  to  zero 
on  the  remaining  three  sides.  These  homogeneous  boundary  conditions  in  the 
x  direction  have  been  accommodated  in  the  present  incremental  formulation 
for  ip  and  w  in  the  second  sweep  of  the  ADI  procedure.  The  point  image 
approach  of  Burggraf1  has  been  used  for  these  derivative  boundary  conditions 
and  was  again  found  to  be  very  satisfactory.  Results  were  obtained  with  30 
step  sizes  in  either  x  and  y  directions,  for  a  value  of  the  Reynolds  Number 
of  100.  They  are  presented  in  figures  1  and  2  as  the  horizontal  velocity 
profile  thru  the  gridpoint  characterized  by  the  maximum  (absolute)  value 

of  the  stream  function,  and  the  contour  plot  of  the  stream  function  itself. 
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In  figure  1  the  results  of  Rubin  et  al.  ’  obtained  by  means  of  a  spline 


approach  using  28  by  28  meshes  in  each  direction  are  shown  for  comparison. 


The  agreement  is  satisfactory  and  provides  further  evidence  of  the  correctness 

of  the  present  approach.  No  quantitative  comparison  is  given  in  figure  2, 

where  the  stream  function  contours  correspond  to  values  of  -0.01,  -0.015, 

-0.020. .. .-0.095  and  agree  reasonably  well  in  values  and  shape  with  previously 
1  22 

published  results  ’ 


C.  Flow  Past  a  NACA  0012  Airfoil 

Flow  past  a  NACA  0012  airfoil  was  finally  considered  in  order  to  test 

the  present  numerical  technique  in  combination  with  the  transformation  of 

Thompson  et  al.^’^.  In  order  to  avoid  the  difficulties  associated  with  a 

sharp  trailing  edge,  this  has  been  smoothed  by  means  of  a  circular  arc, 

see  figures  3  and  4  which  also  provide  the  transformed  coordinate  lines  in 

the  physical  plane.  The  coordinate  transformation,  as  used  in  the  present 

21 

study,  has  been  kindly  provided  by  Captain  H.  A.  Hegna  who  has  optimized 
the  spacing  of  the  (n  =  constant)  coordinate  lines  near  the  body  sur¬ 
face  for  turbulent  flow  at  a  value  of  the  Reynolds  number  of  about  10^ . 

No  attempt  was  made  in  the  present  study  at  optimizing  the  above  mentioned 
coordinate  spacing  for  the  present  laminar  flow  calculations.  Figure  3 
clearly  shows  the  outer  boundary  of  the  integration  domain  in  the  physical 
plane,  that  is,  a  circle  of  center  in  the  origin  and  radius  equal  to  10. 
The  body  and  the  coordinate  lines  immediately  around  it  are  instead  very 
poorly  resolved  by  the  large  scale  of  the  computer  plot.  Figure  4  shows  a 
blow-up  of  the  airfoil  whose  nondiraensional  chord  length  equals  1  and  of 
the  coordinate  lines  immediately  surrounding  it.  However,  even  such  a 
fairly  large  scale  is  not  sufficient  to  clearly  show  the  spacing  of  the 
n  =  constant  coordinate  lines  immediately  surrounding  the  airfoil. 


r 


Solutions  have  been  obtained  for  flow  at  zero  angle  of  attack  for 

2  4 

two  values  of  the  Reynolds  number,  namely  Re  =  10  and  Re  =  10  .  The 
corresponding  velocity  vector  plots  are  given  in  figures  5  and  6,  respec¬ 
tively.  Figure  5  shows  a  typical  attached  viscous  flow  configuration 
with  a  clearly  visible  boundary  layer  near  the  surface  of  the  body.  In 
figure  6,  the  higher  Reynolds  number  is  clearly  seen  to  produce  a  much 
thinner  boundary  layer.  A  blow-up  of  the  velocity  distribution,  very  close 
to  the  body  surface,  is  given  in  figures  7  and  8  for  the  same  two  flow  con¬ 
figurations.  These  figures  clearly  show  that  despite  a  coordinate  spacing 
independent  of  the  Reynolds  number,  the  numerical  solution  has  been  able 
to  capture  the  shrinking  of  the  boundary  layer  thickness  with  increasing 
Reynolds  number  of  the  flow.  It  is  obvious,  however,  that  an  optimization 

of  the  coordinate  spacing  is  warranted  in  order  to  obtain  highly  accurate 

4  2 

solutions  at  Reynolds  number  values  of  10  or  higher.  A  Re  =  10  flow  at 
an  angle  of  attack  of  0.1  has  also  been  computed  and  the  results  are  given 
in  figure  9  again  as  velocity  vectors.  Separated  flows  were  not  attempted 
due  to  the  fact  that  all  convective  terms  in  the  governing  equations  are 
represented  by  central  finite  differences.  However,  first  order  accurate 
windward  finite  difference  representations  for  such  terms  can  be  easily 
accommodated  in  the  present  algorithm.  The  present  approach  is  computa¬ 
tionally  very  fast  insofar  as  the  solution  proceeds  thru  100  time  steps 
within  about  2  CPU  minutes  of  CDC  Cyber  175  for  the  present  calculations 

g 

employing  a  grid  of  70  by  44  points.  The  method  of  Thompson  et  al 
requires  a  comparable  amount  of  computer  time  for  advancing  the  solution 
of  a  single  time  step.  However,  the  convergence  rate  was  found  to  be  lower 
than  anticipated:  whereas,  for  the  driven  cavity  problem  150  time  steps 
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(At  =  .1)  were  sufficient  for  a  satisfactory  convergence  ([(wn+'*'  -  ujn) 

/wn]  <  10  "*) ,  for  the  flow  past  the  NACA  0012  airfoil  the  relative 

average 

-4 

error  in  the  solution  was  still  of  the  order  of  10  ,  after  200  time  cycles. 

Further,  in  order  to  obtain  convergence,  the  solution  had  to  be  started 

-4 

with  a  very  small  step  size  (At  =10  ) ,  which  was  then  increased  at  every 

_2 

iteration  by  a  factor  of  1.1,  until  a  value  of  10  was  reached,  and  was 
then  kept  constant  at  that  value,  in  order  to  avoid  divergence.  All  results 
have  been  obtained  with  the  backward-in-time  approach.  The  program  using  a 
Crank  Nicolson  averaging  was  not  found  to  converge  for  comparable  values  of 


the  time  step. 


SECTION  V 


CONCLUSIONS  AND  RECOMMENDATIONS 

An  algorithm  for  solving  the  vorticity-stream  function  Navier-Stok.es 

equations  for  steady  laminar  flow  past  an  arbitrary  airfoil  has  been 

developed.  The  governing  equations  are  written  in  a  system  of  body  oriented 

coordinates^* ^  and  solved  by  means  of  the  ADI  procedure  of  Douglas  and 
17 

Gunn  .  The  present  approach  has  computed  the  flow  field  past  a  HACA  0012 

airfoil  successfully,  and  has  shown  to  be  cost-competitive  with  other 

approaches  available  in  the  technical  literature.  Further,  it  could  be 

easily  modified  to  extend  its  capability  to  unsteady  flow  computations  by 

relaxing  the  stream  function  equation  at  every  time  step  by  any  suitable 

numerical  technique  (line  SOR,  ADI,  Direct  solvers).  However,  the  present 

approach  needs  improvements  with  respect  to  the  convergence  rate  and  the 

present  inability  to  compute  separated  flows.  In  this  respect,  the  effect 

A  13 

of  windward  differencing  and  variable  time  steps  ’  is  certainly  worth 

investigating.  Finally,  it  is  the  author's  belief  that  further,  dramatic 

improvements  could  be  obtained  by  incorporating  the  very  promising  multi- 
23 

grid  idea  of  Brandt 
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COORDINATE  LINES 


Figure  3.  NACA  0012  Airfoil;  Body-Oriented  Coordinates  (Far  Field) 
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Figure  8.  NACA  0012  Airfoil;  Re=10,000  Boundary  Layer  Profiles 


APPENDIX 


AN  ALGORITHM  FOR  THE  SOLUTION  OF  PERIODIC  SYSTEMS  OF 


TWO  COUPLED  TRIDIAGONAL  EQUATIONS 


Algorithm 

18 

The  present  Appendix  generalizes  the  algorithm  of  Ahlberg  et  al.  , 
for  solving  tridiagonal  periodic  systems,  to  the  case  of  two-by-two  block 
tridiagonal  systems. 

The  most  general  tridiagonal  system  for  21  coupled  equations  in  21 
unknowns  k^,  F^,  i  =  1,  2  .  .  .1,  with  periodic  boundary  conditions  is 
given  as: 


a^  kj  +  kx  +  cl1  k2  +  dl1  Fz  +  e^  F1  +  fl1  ?2  =  h^ 


a21  kx  +  b21  kx  +  c21  k2  +  d2x  Fj  +  e21  F1  +  f2x  F2  =  h21 


3li  ki-l  +  bli  ki  +  cli  kl+l  +  dll  Fi-1  +  eli  Fi  +  fli  Fi+1  ”  hli  (2a> 

a2l  ki-l  *  b2i  ki  +  a2i  ki+l  +  d2i  Fi-1  +  a2i  Fi  +  £2i  Fi+1  *  b2i  <2b> 


alI  kI-l  *  blI  kI  +  clI  kl  +  dlI  FI-1  +  alI  FI  +  £1I  F1  '  bF! 
a2I  kI-l  +  b2I  kI  +  c2I  kl  +  d2I  FI-1  +  a2I  FI  +  £2I  F1  '  h2f 


Let  us  assume  the  following  recursion  relations  for  the  unknowns,  k^,  F^ : 
ki  =  rli  ki+l  +  sli  Fi+1  +  ai  +  uli  kl  +  vli  Fi  • 


Fi  r2i  ki+l  +  s2i  Fi+1  +  t2i  +  u2i  kI  +  v2i  Fl 


Equations  (4a)  and  (4b)  are  valid  for  any  value  of  i.  They  can  be  used 


therefore,  to  eliminate  k^_1  and  from  equations  (2a, b) ,  which  become: 

bli  ki  +  cli  ki+l  +  eli  Fi  +  fli  Fi+1  +  mli  kI  +  nli  FI  =  hli  ’ 

b2i  k1  +  c2  ki+1  +J2i?i  +  f21  F1+1  +  kj  +  n^  Fj  =  h2±  , 


with: 


bli  “  ali  rli_i  +  dli  r2i-l  +  bli  * 


cli  »  cl1  , 


eli  =  ali  sli-l  +  dli  s2i-l  +  eli  ’ 


fli  =  fli  , 


(5a) 

(5b) 

(6a) 

(7a) 

(8a) 

(9a) 


and : 


mi 

i 

=  al . 
l 

uli-i 

+  dli 

u2i-l 

> 

(10a) 

nTi 

=  a  l  . 
i 

vli-i 

+  dlt 

v2i-l 

* 

(11a) 

bli 

*  hi. 
i 

-  al  . 

l 

lll-l 

-  dli 

t2i-l  ’ 

(12a) 

hi . 

1 

=  a  J  . 
l 

rli-l 

+  d  2  . 

1 

r2l-l 

+  b2 . 
i  , 

(6b) 

1  -  , 

1 

-  .  ) 

-  t  -  , 

t 

> 

(7b) 

7l  . 

1 

=  ul. 

l 

Sli-1 

+  d2l 

s2i-l 

+  e2i  ’ 

(8b) 

m2,  =  a2,  ul,  ,  +  d2,  u2 


'i 


l  1-1  i  i-1  ’ 


n2.  =  a2,  vl,  .  +  d2.  v2 


i  i-1  i  i-1  ’ 


h?i  =  h2 i  -  a2t  tl^  -  d2i  t2t  ^ 


(9b) 

(10b) 

(lib) 

(12b) 
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m. 


The  rli  thru  tli  coefficients  can  now  be  determined  by  multiplying  equations 
(5a)  and  (5b)  by  7J±  and  e^  respectively,  subtracting  them  to  eliminate 
and  solving  for  k^,  to  give: 


rli  - 

(eli 

c2i 

'^i 

*i 

-b2i 

eli). 

sli  “ 

(iIi 

«i 

-  72  ± 

fl^/Cbl^ 

^i 

-bli 

elj^)  , 

uli 

(eli 

m2i 

-^i 

/  (bii 

^i 

““i 

ili). 

Vli  = 

(iri 

^i 

-  12  i 

nl±)/ibi± 

^i 

"“i 

eT^, 

tl.  = 
1 

(^ 

“i 

-ili 

h2_1 )  /  (bl^ 

7*± 

-  “i 

iTp. 

The  r2^  thru  t2^  coefficients  are  similarly  obtained  by  elminating  and 
solving  for  to  give: 

r2t  =  (bl.  72±  -  b2±  7ri)/(Ili  b2.  -  l2±  bl±) , 

s2i  =  (bli  £2i  -  b2L  fTi)/(eTi  b2±  -  i2±  bl±) , 

u2^  =  (bl^  m2^  -  b2^  ml^)/(elj,  b2^  -  e2^  bl^) , 

v2t  =  (bl1  ^2.  -  b2.  ^I.)/(ili  b2i  -  72  ±  bl1), 

t2i  =  (^i  ^i  "  ^i  **V/(®*i  ^i>* 


(13a 

(14a 

(15a 

(16a 

(17a 


(13b 

(14b 

(15b 

(16b 

(17b 


The  rlj  thru  t2^  coefficients  are  evaluated  by  means  of  equations  (la,b)  in 
the  same  way  and  all  the  rl^  thru  t2^  (i=2...,I-l)  can  then  be  evaluated  by 
means  of  equations  (13a-17b) .  In  order  to  evaluate  the  F^  and  k^  unknowns 
(i.e.  to  solve  the  original  system  of  equations)  let  us  now  take: 

ki  =  wuli  kI  +  wvli  Fi  +  wtli»  (18a 

F.^  =  wu2i  kj  +  wv2^  Fj  +  wt2if 

I 


(18b 


with  the  wulj  thru  wt2^  coefficients  obviously  given  as: 

wulj  =  wv2j  =  1,  wvlj  =  wtlj  =  wu2j  =  wt2j  =  0.  (19a, b) 

From  equations  (18a, b)  and  (4a, b)  it  is  easy  to  verify  that: 


wul^  =*  rl^  wuli+^  +  sl^  wu2^+^  +  ul^. 


wvl^  =  rl^  wl^  +  sl^  wv2i+^  +  vl^. 


wtl1  =  rli  «tli+1  +  sl±  wt21+i  +  t;li* 


wu2.  =  r2,  wul.,,  +  s2  wu2  .  +  u2 . , 
l  i  i+l  i  l+l  l 


wv2i  =  r21  wvli+1  +  s2.  wv  1+]_  +  v2., 


wt2.  =  r2.  wtl. , ,  +  s2  wt2  +  t2  . 
i  l  i+l  i  i+l  i 


(20a) 

(21a) 

(22a) 

(20b) 

(21b) 

(22b) 


Equations  (20a-22b)  together  with  "boundary  conditions"  (19a, b)  allow  the 
evaluations  of  the  wu^  thru  wt2^  coefficients  (i  =  1-1,  1-2,..., 1).  These 
will  finally  provide  the  solution  vectors  f^,  k^  if  and  k^  can  be  some¬ 
how  determined.  This  is  easily  accomplished  by  eliminating  the  F^,  k^, 

Fj_^,  kj_^  unknowns  in  equations  (3a, b)  by  means  of  the  appropriate  recursion 
formulas  (18a, b)  and  by  solving  the  resulting  system  of  2  equations  in  2 
unknowns,  k^.,  F^  by  means  of  the  Kramer's  rule. 


Fortran  Implementation 


The  listing  of  a  Fortran  subroutine  implementing  the  present  algorithm 
is  attached  for  convenience.  Note  that  the  possibility  of  using  the  same  arrays 
for  the  recusion  coefficients  ul^,  wul^  .  .  .  .  ,  t2^,  wt2^,  has  been  exploited. 
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subroutine  *er:n9(jenoi 


0)0100 
0  0  QUO 

•  THE  FOLLOHING  SUBROUTINE  IN/ERTES  A  CQUPLEO  SET  OF  T  NO  ••  030120 

•  TRIOIAGONAL  EQUATIONS  WITH  PERIODIC  JOUNOART  CONDITIONS.  •»  010130 

•  THE  COEFFICIENTS  OF  THE  DIFFERENCE  EQUATIONS  ARE  A1  THRU  ••  OJOIOO 

•  HI  AND  A*  THRU  HI.  THE  SOLUTION  9ECTJRS  ARE  RETURNED  TO  ••  010150 

•  THE  MAIN  PROGRAN  IT  MEANS  OF  THE  SECTORS  HI  AND  HZ.  «♦••**  jOGlbO 

. . . . .  0  3  0170 

COMMON  AltiSJ)  ,  Jl  (40)  ,D1(80I,D1  1601  ,E1  (80)  ,F11301  .HllOO)  ,  0001(0 

1A2 (831  ,  82<8CI,:2<301,Q2l(i(j>,E2<eJI,F2(80l,H2(501  00  0190 

OIHENSION  RK4U  ,Sl<tCI  ,  TU40  )  ,U1<  10*  »J1<  801 .  0JC20C 

1R2I03I  , S2I631 ,  12  (90  ,  JZ(ST) ,  02  <801  00  0210 

JMUJENQ-1  C10220 

OENCi‘1 ./(91«1)*E2I1)-02C1»»E1(1)»  000230 

RKIDOENON*  <E1<  1)  *C2<  1)  -E2<  It  *C1I  II  >  0  ]  0  2*0 

Sllll *DEN0M*<E1<  1» •F2<1I-E2(D*F1(1I1  03 t 250 

T1  ( 11  =OENOM*  <S2(11*H1<11-E1<U*H2(111  0  3  0260 

Ul(l)  *OENOH*  CElt  1!*A?<  II  -E2I1  >*A1C  ID  00  0270 

7K1MJENJH*  <E1  (  11*D2(  11  -E2(U»C1<  1)  1  3  3  0280 

R2  I 11  =  DENO‘(*  1 1.  ll>»CUU-eUU*C2<lD  000290 

S2<ll=OENOM*  <32(11 ‘Fll  ll-ei<ll»f2(lll  00030  0 

T2 <  11  «OENOM*  UK  11»H2(  11  -  62(11  •HUH  1  00  0310 

U2(1»*DEnOH*(3E(U»A1<  11-61(1 1  ‘32(111  0JC320 

92<11=0EN0H‘  (Id  (11  •01(11-81(11  *02(  II  I  QO  0330 

DO  6  J=2,JM1  00  0340 

91N*41(J1»R1(  J-1D01  <J1*R2(J-1I*BKJ1  0  3G350 

C1N*C1 ( J1  0  j 0360 

E1N*A1  ( Jl*Sl  (J-1H01  ( J  1  •  S2  ( J  - 1 1  *E1(J1  0  3  0370 

FINiFKJI  000360 

AM1N=A1(J1»J1(J-11  ♦01<  J1  »U2<  J-U  0  3  0390 

AN1N=A1 (JI*J1(J-11*DKJI*92(J-U  000400 

H1N-M1 ( J1 -Al(Jl  •T1<J-11*D1(JI  *72(8-11  3  3  0410 

B2N*A2(JI»Rlt  J-U*02(J1*R2IJ-1I  *82(31  C3042D 

C2N*C2(J1  0  3  C  430 

E2N*A2<  JI»St<  J-U*02(JI  *S2  (J-U  *£2  (  Jl  0  3  0440 

F2n*F2 ( J1  030450 

AH2N=A2<JI»J1<J-1I»02<  J1  »U2(  J-1)  0  3  0460 

AN2N=A2  (  Jl  »;i  (  J-1I*D2(  JI  *J2<  J-U  3  3  0470 

M2N*H2  <JI-A2(JI*T1  (J-1I-02(JI»T2(J-U  03  0460 

DEN*1 ./ <81N*t2N-92N*ElNl  030490 

R1  ( Jl  »0EN»(E1N*C2N-E2N»C1NI  0  0  0500 

Sl< J1*0EN*(E1S‘F2N-E2N*F1N1  010510 

Tl(  Jl  *OEN*  (E2.N*N1N-E1N*H2N1  0  0  05  20 

U1  (Jl  IQEN*  (E1N*  AF2N-E2N»AH1N)  0  3  0530 

91<JlxO£N*(ElN*AN2N-t2N*ANlNI  030540 

R2(  J1*-UEN*(31U»C2N-E2N*;1N|  0  3  05  50 

S2<  JM-OEN*  <8i:i*  F2N-r?2N»FlNI  00C560 

T2 (J1»-0£n»<Q2n»hin-rin»M2ni  030570 

U2<  Jl  «-0EN*(31N*AM2N-R2N*Am1NI  030580 

W2< Jl s-OEN*(31N*AN2N-32N*  AN1N1  0  3  0590 

6  CONTINUE  030600 

UKJENOUl.  000610 

71 ( JEND1 s0 ■  030620 

T1(JENOI*0.  030630 

U2IJENDU0.  000640 

V2 ( JENO I *1 •  030650 

T2(JEN01>0.  030660 

DO  7  J*l,JHl  03067C 

<■ JENO- J  030600 

U1(K1*U1<K1*RI(<I»U1(<*1I *SK<I*U2<«*11  0  3  0690 

71(K16J1(KI*R1(<I»/1(<*11*S1(K1*720(*11  03  0700 

T1(K1 " T 1 <KI *R1 ( <1 • T1 <5*1 1 *S1 INI *T2 ( K  ♦  1 1  030710 

U2 (<1 *U2 (51 *R2i<l*UI (4*1 1 *S2t<l *U2 <K*it  030720 

V2<KI*02(K1*R2<<I*91<<»11 *S2(KI»92(K*1I  000730 

T2(K1  «T2(K1*R2(K1»T1  (K*u  *S2(KI»T2«*11  030740 

7  CONTINUE  000750 

N*  J  END  030760 

H*  JHi  000770 

AlFIXI  (N1*U1  (HI  *91(NI  *C  1  (Nl  *Ul  (It  *01  (NI*U2(NI  *F1  (N1*U2(  1 1  0-'  0  700 

BETl*  A1  <Nt*9l  <  II  *E1  (Nl  *C1  (Nl  *91  <11  ♦OKNI  •V2CM  *FltM*92(  11  0  3  0790 

GAMAHHl (N) -Ai (NI»Tl (*l -l1(ni >11(1 1-01 INI •T!(1I-F1(NI  »T2  <  11  030600 

ALF2*A2(NI»U1(H1  *92 (til  *C2  (Nl  »’J  l  ( 11  *02  (Nl  •U2IM1  *F2  (N)»  J2<  11  000810 

8ET2*A2  (Nl* Oil Nl *E2( Nl *32 (Nl ’ 01 (11 *02(Nl»»2(il *F2 (Nl >92(1  1  0  3  0820 

GAMA2»h2IN1-A2(NI»T1 I*1-C2(N1 • Tltll -02 (Nl »T2< N1-F2 (Ml»T2 ( 11  03 0630 

0£N*1./  (ACFl*9rT  Z-ALF?*fU  Til  00  004  0 

Hl<  JENOl  =JEN*<  ..AHA  t»RE  T2-&ANA2*  FET1I  G3C650 

H2 ( JEHU  I •DEN* (GAHA2» A* FI - &AHA 1 • ALF2 1  0  0  00  60 

DO  8  J* 1 i JNl  030670 

HI ( Jl *U1 I J) *H1 ( JTN31 *91 ( Jl *H2 ( JtNOI *T1( Jl  000060 

H2( JI*U2(J1*M1 1 JENOI *92 ( JI*H? ( JEN01 * T 2 ( Jl  0  J  0690 

(  CONTINUE  000900 

RETURN  030910 

END  000920 
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Numerical  Application 

20 

The  present  algorithm  has  been  applied  to  a  Spline  4  numerical 
solution  of  the  following  ordinary  differential  equation: 

F" (x)  +  F (x)  =  (1  -  4n2)  sin  2nx  ,  (23) 

whose  exact  solution,  F(x)  *  sin  2ttx,  is  available  for  comparison.  The 
20 

Spline  4  procedure  ,  applied  to  the  numerical  integration  of  equation  (23) 
in  the  range  0  x  <  1  leads  to  a  system  of  two  coupled  tridiagonal  equations 
(la-3b)  with  the  following  values  of  the  coefficients: 

ali  -  cli  =  1/12;  blt  »  5/6;  dl±  =  fl±  -  0;  el±  =  1; 

hl^  =  (1  -  4ti2)  sin  2tt ( i  -  l)h;  (24a-e) 

a2^  =  c2^  =  h2/6;  b2^  =  2h2/3;  d2^  =  f2^  «  -1;  e2^  =  2;  h2^  =  0.  (25a-e) 

20 

F^  and  are  the  functional  and  second  derivative  Spline  4  values  at 

the  nodal  point  x^tx^  -  (i-l)hj  and  h  is  the  step  size  (h=l/I). 

Numerical  solutions  have  been  obtained  for  four  values  of  I  and  the  cor- 

I 

responding  average  truncation  errors  e(e  =  l  |  F.  -  sin  2n (i-l)h | )/I) ,  are 

1=1  1  4 

given  in  Table  I.  The  errors  are  proportional  to  h  ,  as  they  should  be,  thus 
verifying  the  validity  of  the  proposed  algorithm.  The  solution  corresponding 
to  I  =  160  required  only  114  ms  of  CDC  Cyber  175  computer. 


N 

20 

40 

80 

160 

e 

.17  •  10'4 

.11  •  io~5 

.69  •  10'7 

00 

1 

o 

r—i 

• 

m 

TABLE  I 
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